
%This is the code to bootstrap the standard errors. --> only works if NCS
%is constant over NP. Meaning each Agent faces the same amount of choice
%situations. 
MeanEst=paramhat;
mnhold=zeros((NF+NV+NV),NReps);
XMAT_Original=XMAT;
PID_Original=XMAT(:,1);
CSID_Original=XMAT(:,2);

clear i r
for xi=1:NReps
    btsample=randi(NP,[NP,1]);
    XMAT=[];
        for r=1:NP
              thisp=btsample(r,1);
              thisx=XMAT_Original(PID_Original==thisp,:);
              thisx(:,1)=r.*ones(size(thisx,1),1);
              XMAT=[XMAT ; thisx];  
        end
        %Make sure that choice situation and ID are identical and going
        %from 1 to NS: 
        %XMAT(:,2)=CSID_Original;
        XMAT(:,2)=XMAT(:,1);
   create_arrays
   create_draws
 tic;
options=optimset('LargeScale','off','Display','iter','GradObj','on',...
    'MaxFunEvals',10000,'MaxIter',MAXITERS,'TolX',PARAMTOL,'TolFun',LLTOL,'DerivativeCheck','off');
[paramhat,fval,exitflag,output,grad,hessian]=fminunc(@loglik,param,options);

disp(' ');
disp(['Estimation took ' num2str(toc./60) ' minutes.']);

    fvalhold(xi) = fval;
    paramhold(:,xi)=paramhat;
    mnhold(:,xi)=paramhat ;   
end
MeanSE=std(mnhold,0,2);



%Segment parameters and account for unestimated zero mean in distribution 5
if NF>0
  fhat=MeanEst(1:NF,1);
  fsd=MeanSE(1:NF,1);
end

if NV>0
  if sum(IDV(:,2) == 5) >0;
     bhat=zeros(NV,1);
     bsd=zeros(NV,1);
     bhat(IDV(:,2) ~= 5,1)=MeanEst(NF+1:NF+sum(IDV(:,2) ~= 5),1);
     bsd(IDV(:,2) ~= 5,1)=MeanSE(NF+1:NF+sum(IDV(:,2) ~= 5),1);
     what=MeanEst(NF+sum(IDV(:,2) ~= 5)+1:end,1);
     wsd=MeanSe(NF+sum(IDV(:,2) ~= 5)+1:end,1);
  else;
     bhat=MeanEst(NF+1:NF+NV,1);
     bsd=MeanSE(NF+1:NF+NV,1);
     what=MeanEst(NF+NV+1:NF+NV+NV,1);
     wsd=MeanSE(NF+NV+1:NF+NV+NV,1);
  end;
end


disp('RESULTS');
disp(' ');
disp(' ');

if NF>0
disp('FIXED COEFFICIENTS');
disp(' ');
disp('                      F      ');
disp('              ------------------ ');
disp('                Est         SE ');
for r=1:length(NAMESF);
    fprintf('%-10s %10.4f %10.4f\n', NAMESF{r,1}, [fhat(r,1) fsd(r,1)]);
end
disp(' ');
end

if NV>0
disp('RANDOM COEFFICIENTS');

disp(' ');
disp('                      B                      W');
disp('              ------------------   -----------------------');
disp('                 Est     SE            Est         SE');
for r=1:length(NAMES);
    fprintf('%-10s %10.4f %10.4f %10.4f %10.4f\n', NAMES{r,1}, [bhat(r,1) bsd(r,1) what(r,1) wsd(r,1)]);
end

end



